
## @knitr MODELS_SETUP_PH

rm(list=ls())

#library(foreign)
#library(stringr)
#library(xlsx)

# from \url{https://www.r-bloggers.com/identifying-the-os-from-r/}
get_os <- function(){
  sysinf <- Sys.info()
  if (!is.null(sysinf)){
    os <- sysinf['sysname']
    if (os == 'Darwin')
      os <- "osx"
  } else { ## mystery machine
    os <- .Platform$OS.type
    if (grepl("^darwin", R.version$os))
      os <- "osx"
    if (grepl("linux-gnu", R.version$os))
      os <- "linux"
  }
  tolower(os)
}

get_machine  <- function(){
    sysinf  <- Sys.info()
    if (!is.null(sysinf)){
        my.machine  <- sysinf['nodename']
    }
    tolower(my.machine)
}


if(get_os()=="windows" & get_machine()=="mompracen") {
    path <- setwd('C:/Users/giaco/Documents/MEGAsync/')
} else if (get_os()=="windows" & get_machine()!="mompracen") {
    path <- setwd('C:/Users/grieco/Dropbox/')
} else if (get_machine()=="aramis") {
    path <- setwd('/home/giacomo68/Dropbox/')
} else {
    path <- setwd('/home/giacomo/Dropbox/')
}


library(survival)
library(xtable)
library(ggplot2)
library(english)
library(mice)
library(mitools)
library(sandwich)
library(lmtest)
library(miceadds)


Data <- read.csv('gc-jg-project/kampala/data/data-kampala.csv',
                 stringsAsFactors=FALSE)


Data$country_id <- as.numeric(factor(Data$country_name))

m <- 10

imp <- mice(Data[,c("kampala",
                 "sidea.1991",
                 "sideb.1991",
                 "my.mil.cap.lag",
                 "v2x_polyarchy.lag",
                 "v2x_libdem.lag",     
                 "v2x_partipdem.lag",
                 "v2x_delibdem.lag", 
                 "v2x_egaldem.lag",
                 "gdp.cap.lag",
                 "pop.lag",
                 "common.law",
                 "year",
                 "rule.law.lag",
                 "aid.recipient",
                 "aid.recipient.q",
                 "All_NGO.lag",
                 "country_id")], m = m, maxit = 10, seed = 1510)

models <- list()
model_results <- list()
merged_imputations <- list()
m1.model_results <- list()
m2.model_results <- list()
m3.model_results <- list()
m4.model_results <- list()
m5.model_results <- list()
m6.model_results <- list()

fit.m1 <- list()
fit.m2 <- list()
fit.m3 <- list()
fit.m4 <- list()
fit.m5 <- list()
fit.m6 <- list()


for (i in 1:m) {
    complete_data <- complete(imp, i)

    complete_data$.imp <- i
    
    merged_data <- merge(complete_data, 
                         Data[,c("country_id","year","country_name",
                                 "kampala.t0","kampala.t1",
                                 "sidea.1945",
                                 "sideb.1945",
                                 "sidea.force.1945",
                                 "sideb.force.1945",
                                 "sidea.no.force.1945",
                                 "sideb.no.force.1945",
                                 "sidea.force.1991",
                                 "sideb.force.1991",
                                 "sidea.no.force.1991",
                                 "sideb.no.force.1991"
                                 )], 
                        by = c("country_id", "year"),  # Note: quoted variable names
                        all.x = TRUE)  # Keep all rows from complete_data
    
    merged_imputations[[i]] <- merged_data

    current_data <- merged_imputations[[i]]

    m1 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
              sidea.1991
            + sideb.1991
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_id)
           ,data=current_data)

     m1.vcov_cluster <- vcovCL(m1, cluster = current_data$country_id)

     m1.model_results[[i]] <- list(
        coefficients = coef(m1),
        variance = m1.vcov_cluster,
        n = nrow(current_data)
       )
    
    fit.m1[[i]] <- cox.zph(m1, transform='log')
    
    m2 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                sidea.1945
            + sideb.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m2.vcov_cluster <- vcovCL(m2, cluster = current_data$country_id)

    m2.model_results[[i]] <- list(
        coefficients = coef(m2),
        variance = m2.vcov_cluster,
        n = nrow(current_data)
       )
    
    fit.m2[[i]] <- cox.zph(m2, transform='log')

    m3 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.force.1991
                + sideb.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)

    m3.vcov_cluster <- vcovCL(m3, cluster = current_data$country_id)

    m3.model_results[[i]] <- list(
        coefficients = coef(m3),
        variance = m3.vcov_cluster,
        n = nrow(current_data)
       )
    
    fit.m3[[i]] <- cox.zph(m3, transform='log')

    m4 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.force.1945
            + sideb.force.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m4.vcov_cluster <- vcovCL(m4, cluster = current_data$country_id)

    m4.model_results[[i]] <- list(
        coefficients = coef(m4),
        variance = m4.vcov_cluster,
        n = nrow(current_data)
       )
    
    fit.m4[[i]] <- cox.zph(m4, transform='log')

    m5 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.no.force.1991
                + sideb.no.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m5.vcov_cluster <- vcovCL(m5, cluster = current_data$country_id)

     m5.model_results[[i]] <- list(
        coefficients = coef(m5),
        variance = m5.vcov_cluster,
        n = nrow(current_data)
       )

    fit.m5[[i]] <- cox.zph(m5, transform='log')

    m6 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.no.force.1945
                + sideb.no.force.1945
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m6.vcov_cluster <- vcovCL(m6, cluster = current_data$country_id)

    m6.model_results[[i]] <- list(
        coefficients = coef(m6),
        variance = m6.vcov_cluster,
        n = nrow(current_data)
       )

    fit.m6[[i]] <- cox.zph(m6, transform='log')
    
}


m1.chisq_values <- sapply(fit.m1, function(x) x$table[, "chisq"])
m2.chisq_values <- sapply(fit.m2, function(x) x$table[, "chisq"])
m3.chisq_values <- sapply(fit.m3, function(x) x$table[, "chisq"])
m4.chisq_values <- sapply(fit.m4, function(x) x$table[, "chisq"])
m5.chisq_values <- sapply(fit.m5, function(x) x$table[, "chisq"])
m6.chisq_values <- sapply(fit.m6, function(x) x$table[, "chisq"])

m1.v1 <- micombine.chisquare(m1.chisq_values[1,], 1, display=FALSE, version=1) #sidea
m1.v2 <- micombine.chisquare(m1.chisq_values[2,], 1, display=FALSE, version=1) #sideb
m1.v3 <- micombine.chisquare(m1.chisq_values[3,], 1, display=FALSE, version=1) #milex
m1.v4 <- micombine.chisquare(m1.chisq_values[4,], 1, display=FALSE, version=1) #polyarchy
m1.v5 <- micombine.chisquare(m1.chisq_values[5,], 1, display=FALSE, version=1) #gdp
m1.v6 <- micombine.chisquare(m1.chisq_values[6,], 1, display=FALSE, version=1) #pop
m1.v7 <- micombine.chisquare(m1.chisq_values[7,], 1, display=FALSE, version=1) #common_law
m1.g1 <- micombine.chisquare(m1.chisq_values[8,], 7, display=FALSE, version=1) #global

m2.v1 <- micombine.chisquare(m2.chisq_values[1,], 1, display=FALSE, version=1) #sidea
m2.v2 <- micombine.chisquare(m2.chisq_values[2,], 1, display=FALSE, version=1) #sideb
m2.v3 <- micombine.chisquare(m2.chisq_values[3,], 1, display=FALSE, version=1) #milex
m2.v4 <- micombine.chisquare(m2.chisq_values[4,], 1, display=FALSE, version=1) #polyarchy
m2.v5 <- micombine.chisquare(m2.chisq_values[5,], 1, display=FALSE, version=1) #gdp
m2.v6 <- micombine.chisquare(m2.chisq_values[6,], 1, display=FALSE, version=1) #pop
m2.v7 <- micombine.chisquare(m2.chisq_values[7,], 1, display=FALSE, version=1) #common_law
m2.g1 <- micombine.chisquare(m2.chisq_values[8,], 7, display=FALSE, version=1) #global

m3.v1 <- micombine.chisquare(m3.chisq_values[1,], 1, display=FALSE, version=1) #sidea
m3.v2 <- micombine.chisquare(m3.chisq_values[2,], 1, display=FALSE, version=1) #sideb
m3.v3 <- micombine.chisquare(m3.chisq_values[3,], 1, display=FALSE, version=1) #milex
m3.v4 <- micombine.chisquare(m3.chisq_values[4,], 1, display=FALSE, version=1) #polyarchy
m3.v5 <- micombine.chisquare(m3.chisq_values[5,], 1, display=FALSE, version=1) #gdp
m3.v6 <- micombine.chisquare(m3.chisq_values[6,], 1, display=FALSE, version=1) #pop
m3.v7 <- micombine.chisquare(m3.chisq_values[7,], 1, display=FALSE, version=1) #common_law
m3.g1 <- micombine.chisquare(m3.chisq_values[8,], 7, display=FALSE, version=1) #global

m4.v1 <- micombine.chisquare(m4.chisq_values[1,], 1, display=FALSE, version=1) #sidea
m4.v2 <- micombine.chisquare(m4.chisq_values[2,], 1, display=FALSE, version=1) #sideb
m4.v3 <- micombine.chisquare(m4.chisq_values[3,], 1, display=FALSE, version=1) #milex
m4.v4 <- micombine.chisquare(m4.chisq_values[4,], 1, display=FALSE, version=1) #polyarchy
m4.v5 <- micombine.chisquare(m4.chisq_values[5,], 1, display=FALSE, version=1) #gdp
m4.v6 <- micombine.chisquare(m4.chisq_values[6,], 1, display=FALSE, version=1) #pop
m4.v7 <- micombine.chisquare(m4.chisq_values[7,], 1, display=FALSE, version=1) #common_law
m4.g1 <- micombine.chisquare(m4.chisq_values[8,], 7, display=FALSE, version=1) #global

m5.v1 <- micombine.chisquare(m5.chisq_values[1,], 1, display=FALSE, version=1) #sidea
m5.v2 <- micombine.chisquare(m5.chisq_values[2,], 1, display=FALSE, version=1) #sideb
m5.v3 <- micombine.chisquare(m5.chisq_values[3,], 1, display=FALSE, version=1) #milex
m5.v4 <- micombine.chisquare(m5.chisq_values[4,], 1, display=FALSE, version=1) #polyarchy
m5.v5 <- micombine.chisquare(m5.chisq_values[5,], 1, display=FALSE, version=1) #gdp
m5.v6 <- micombine.chisquare(m5.chisq_values[6,], 1, display=FALSE, version=1) #pop
m5.v7 <- micombine.chisquare(m5.chisq_values[7,], 1, display=FALSE, version=1) #common_law
m5.g1 <- micombine.chisquare(m5.chisq_values[8,], 7, display=FALSE, version=1) #global

m6.v1 <- micombine.chisquare(m6.chisq_values[1,], 1, display=FALSE, version=1) #sidea
m6.v2 <- micombine.chisquare(m6.chisq_values[2,], 1, display=FALSE, version=1) #sideb
m6.v3 <- micombine.chisquare(m6.chisq_values[3,], 1, display=FALSE, version=1) #milex
m6.v4 <- micombine.chisquare(m6.chisq_values[4,], 1, display=FALSE, version=1) #polyarchy
m6.v5 <- micombine.chisquare(m6.chisq_values[5,], 1, display=FALSE, version=1) #gdp
m6.v6 <- micombine.chisquare(m6.chisq_values[6,], 1, display=FALSE, version=1) #pop
m6.v7 <- micombine.chisquare(m6.chisq_values[7,], 1, display=FALSE, version=1) #common_law
m6.g1 <- micombine.chisquare(m6.chisq_values[8,], 7, display=FALSE, version=1) #global


m1.chi2 <- rbind(m1.v1,
                 m1.v2,
                 m1.v3,
                 m1.v4,
                 m1.v5,
                 m1.v6,
                 m1.v7,
                 m1.g1)

m2.chi2 <- rbind(m1.v1,
                 m2.v2,
                 m2.v3,
                 m2.v4,
                 m2.v5,
                 m2.v6,
                 m2.v7,
                 m2.g1)

m3.chi2 <- rbind(m3.v1,
                 m3.v2,
                 m3.v3,
                 m3.v4,
                 m3.v5,
                 m3.v6,
                 m3.v7,
                 m3.g1)

m4.chi2 <- rbind(m4.v1,
                 m4.v2,
                 m4.v3,
                 m4.v4,
                 m4.v5,
                 m4.v6,
                 m4.v7,
                 m4.g1)

m5.chi2 <- rbind(m5.v1,
                 m5.v2,
                 m5.v3,
                 m5.v4,
                 m5.v5,
                 m5.v6,
                 m5.v7,
                 m5.g1)

m6.chi2 <- rbind(m6.v1,
                 m6.v2,
                 m6.v3,
                 m6.v4,
                 m6.v5,
                 m6.v6,
                 m6.v7,
                 m6.g1)




ph.test1 <- round(rbind(cbind(m1.chi2[,c(1,3,2)],m3.chi2[,c(1,3,2)],m5.chi2[,c(1,3,2)]),
                        cbind(m2.chi2[,c(1,3,2)],m4.chi2[,c(1,3,2)],m6.chi2[,c(1,3,2)])),3)

rownames(ph.test1) <- rep(c("MID side A",
                            "MID side B",
                            "Milex per cap",
                            "Polyarchy",
                            "GDP cap (log)",
                            "Population (log)",
                            "Common law",
                            "Global test"),2)


ph.test1d <- data.frame(cbind(Variables=rownames(ph.test1),ph.test1))
colnames(ph.test1d) <- c("Variables",rep(c("chisq","df","p"),3))
rownames(ph.test1d) <- NULL


addtorow.ph1 <- list()
addtorow.ph1$pos <- list()
addtorow.ph1$pos[[1]] <- c(-1) # toprule
addtorow.ph1$pos[[2]] <- c(-1) # multicol
addtorow.ph1$pos[[3]] <- c(-1) # midrul
addtorow.ph1$pos[[4]] <- c(8)
addtorow.ph1$pos[[5]] <- c(8)
addtorow.ph1$pos[[6]] <- c(8)
addtorow.ph1$pos[[7]] <- c(nrow(ph.test1d)) # endrule

addtorow.ph1$command <- c('\\toprule\n',
                      "& \\multicolumn{3}{c}{Overall, from 1991}
                       & \\multicolumn{3}{c}{Force, from 1991}
                       & \\multicolumn{3}{c}{No force, from 1991} \\\\\n",
                      '\\midrule\n',
                      '\\midrule\n',
                       "& \\multicolumn{3}{c}{Overall, from 1945}
                       & \\multicolumn{3}{c}{Force, from 1945}
                       & \\multicolumn{3}{c}{No force, from 1945} \\\\\n",
                      '\\midrule\n',
                      '\\bottomrule')


## coup.xtable <- xtable(coup.table1,
##                       caption=paste0("Coups and Attempted Coups in US Partner Nations, ",toString(min.year),"-",toString(max.year)),
##                       label="tab:coups")
## align(coup.xtable) <- "llcccclcccclccc"
## names(coup.xtable) <- NULL
## 


kampala.ph1.plus <- xtable(ph.test1d,
                          caption="Proportional Hazard Tests for Models in Figure~\\ref{fig:simple-regression}",
                          label="tab:ph-tests")
align(kampala.ph1.plus) <- "rlrrrrrrrrr"


print(kampala.ph1.plus,sanitize.colnames.function=identity,
      sanitize.rownames.function=identity,
      sanitize.text.function = function(x) {x},
      add.to.row = addtorow.ph1, include.colnames = TRUE,
      floating = TRUE,
#      floating.environment = "sidewaystable",
      booktabs = TRUE,
      table.placement = "!ht", #"!htbp",
      caption.placement = "top",
      include.rownames=FALSE,
      hline.after=NULL,
      scalebox=.89,
      size="\\normalsize",
     ,file="gc-jg-project/kampala/tables/ph-test1.tex"
      )
